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Detection  of  sub- surface  optical  layers  in  marine  waters  has  important  applications  in  fisheries  management, 
climate  modeling,  and  decision-based  systems  related  to  military  operations.  Concurrent  changes  in  ihe 
magnitude  and  spatial  variability  of  remote  sensing  reflectance  (Rn)  ratios  and  submerged  scattering  layers 
were  investigated  in  coastal  waters  of  the  northern  Gulf  of  Alaska  during  summer  of  2002  based  on  high 
resolution  and  simultaneous  passive  (MicroSAS)  and  active  (Fish  Lidar  Oceanic  Experimental.  FLOE)  optical 
measurements.  Principal  Component  Analysis  revealed  that  the  spatial  variability  of  total  lidar  backscattering 
signal  (S)  between  2.1  and  20  m  depth  was  weakly  associated  with  changes  in  the  inherent  optical 
properties  (lOPs)  of  surface  waters.  Also  based  on  a  250-m  footprint,  the  vertical  attenuation  of  S  was 
inversely  related  to  the  lOPs  (Spearman  Rank  Correlation  up  to  0  43).  Low  (arithmetic  average  and 
standard  deviation)  and  high  (skewness  and  kurtosis)  moments  of  Rn(443)/Rf,(490)  and  Rn(508)/R„(555) 
ratios  were  correlated  with  vertical  changes  in  total  lidar  backscattering  signal  (S)  at  different  locations.  This 
suggests  the  use  of  sub-pixel  ocean  color  statistics  to  infer  the  spatial  distribution  of  sub-surface  scattering 
layers  in  coastal  waters  characterized  by  stratified  conditions,  well  defined  S  layers  (i.e.,  magnitude  of  S 
maximum  comparable  to  near  surface  values),  and  relatively  high  vertically  integrated  phytoplankton 
pigments  in  the  euphotic  zone  (chlorophyll  a  concentration  >150  mg  m  2). 

©  2010  Elsevier  Inc.  All  rights  reserved. 


1.  Introduction 

Remote  sensing  reflectance  above  the  sea  surface  (Rn  (A,  Of), 
where  A  is  the  wavelength)  and  total  lidar  backscattering  signal  (S 
(z))  are  theoretically  connected  in  two  ways.  First,  both  depend  on  the 
volume  backscatter  function  ft  (0),  where  6  is  the  scattering  angle. 
is  directly  proportional  to  backscattering  coefficient  (bb  (A)),  which  is 
an  inherent  optical  property  (IOP)  equal  to  the  integral  of  ft  over  all 
scattering  angles  greater  than  0  5n.  The  amplitude  (A^  (z))  of  the  total 
volume  lidar  backscattering  is  proportional  to  ft  (n).  The  depth 
dependence  of  the  lidar  measurement  is  given  by 

S(z)  =  A<(z)e'-2“Z)  (1) 

where  A^z)  is  the  sum  of  the  volume  backscattering  (ft)  contributions 
due  to  particulate  (e.g.,  phytoplankton,  zooplankton,  fish)  and  dis¬ 
solved  (e.g.,  seawater)  components,  and  a  is  the  vertically  attenuation 
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coefficient  of  lidar  volume  backscattering  and  independent  of 
depth.  Note  that  a  is  equal  to  the  vertically  diffuse  attenuation  of 
downward  irradiance,  Kd,  if  the  lidar  field  of  view  is  relatively  wide 
(Gordon,  1982).  The  17mrad  field  of  view  of  the  lidar  generally 
satisfied  this  condition  (Churnside  et  al.,  1998).  Second,  Rn  is  inversely 
related  to  a.  the  total  absorption  coefficient  (dissolved  and  particulate 
optical  constituents),  an  IOP  with  major  influence  on  Kd,  and  therefore 
responsible  for  lidar  backscattering  attenuation  along  the  vertical 
component 

Current  techniques  to  derive  vertical  shape  of  optical  constituents 
based  on  Rn  have  large  uncertainties  ( — 70%)  or  require  prior 
knowledge  of  depth  weighting  functions  (Piskozub  et  al.,  2008: 
Zanaveld,  1982).  As  stated  by  Gordon  and  Brown  (1975),  the  use  of  R^ 
(A,  Of )  spectra  to  infer  vertical  profiles  of  lOPs  is  impossible  without  a 
pnon  information.  For  example,  even  in  the  simple  case  of  a  single 
sub-surface  peak  of  lOPs  with  a  Gaussian  profile,  not  all  fitting 
parameters  (i.e.,  amplitude,  depth  and  width  of  peak)  could  be 
retrieved  based  solely  on  Rn  (A,  Of  )  measurements. 

The  use  of  Look-Up  Tables  (LUTs)  of  regionally  and  seasonally 
averaged  lOPs  (i.e.,  weighting  averages),  including  dependence  of 
depth,  has  been  one  approach  to  use  surface  observations  of  Rn- 
derived  chlorophyll  a  concentration  (chi  a)  to  obtain  information 
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about  the  vertical  structure  of  phytoplankton  light  absorption 
pigments  in  marine  waters  (Platt  et  al.,  1988;  Sathyendranath  et  al., 
1995,  hereafter  SP  method).  In  this  case,  different  vertical  profiles  of 
lOPs  are  derived  from  ship-based  chi  o  measurements  by  fitting  to 
Gaussian  curves.  This  technique  has  been  a  step  forward  to  elucidate 
large-scale  differences  in  the  vertical  structure  of  lOPs  in  case  I  waters 
(i.e.,  those  where  optical  properties  are  defined  by  phytoplankton  and 
covarying  components).  However,  the  SP  method  does  not  exploit 
other  spatial  statistical  metrics  (e.g..  standard  deviation,  skewness, 
kurtosis)  beyond  the  calculation  of  the  arithmetic  average  of  vertical 
distribution  of  lOPs  and  associated  ocean  color  signatures  derived 
from  remote  sensors. 

The  higher  moments  related  to  “peakedness"  (kurtosis)  and 
“symmetry”  (skewness)  of  the  probability  distribution  have  been 
used  in  the  ocean  color  community  to  characterize  surface  peaks  of 
light-absorbing  particulates  (e.g.,  skewness  and  kurtosis  are  relatively 
high  in  Trichodesmium  blooms)  (Westeberry  &  Siegel.  2006).  Based 
on  an  analysis  of  the  vertical  microscale  distributions  of  fluorescence, 
Mitchell  et  al.  (2008)  concluded  that  kurtosis  and  skewness  of 
phytoplankton  patch  patterns  were  determined  by  multiple  mechan¬ 
isms  (e.g.,  active  migration  of  cells  and  formation  of  organic 
aggregates).  Therefore,  for  a  population  of  pixels,  the  statistical 
distribution  of  Rn  (K  04- )  ratios  can  potentially  provide  information 
about  underwater  perturbations  of  light  fields. 

In  this  study  we  performed  aerial  surveys  using  coincident  passive 
and  active  optical  measurements  in  coastal  waters  of  the  northern 
part  of  the  Gulf  of  Alaska  (NGOA)  in  order  to  answer  the  following 
scientific  questions;  How  does  the  vertical  structure  of  optical 
properties  affect  the  spatial  statistics  of  passive  ocean  color  data 
measured  above  the  ocean  surface?  Is  a  single  lidar  frequency 
sufficient  to  allow  a  ’bridge’  with  above- water  radiance  observations 
collected  by  multi-wavelength  passive  optical  systems? 

2.  Methods 

2.  J.  Study  oreo 

Eastern  shelf  waters  of  Afgonak/Kodiak  Islands  (57.48°-58.04’  N, 
152.91°-!  51.67°  W,  Fig.  1)  during  late  summer  (August  17,  2002)  are 
characterized  by  a  well-developed  pycnocline  (median  vertical  density 
difference  =  0.025  kgm  4  between  0  and  20  m  in  depth)  (Appendix  A, 


Fig.  1.  Study  area  and  aerial  surveys.  Concurrent  ROE  and  MicroSAS  measurements 
(solid  dots),  case  studies  with  HS  (green  dots)  and  IS  (red  dots)  lidar  waveforms 
(Table  3).  Lidar  data  shown  in  Fig.  2  (blue  cross)  1  Afgonak  Island,  2:  Kodiak  Island, 
3:  Marmot  Bay,  4:  Whale  Island.  5:  Spruce  Island.  6.  Chiniak  Bay 


Fig.  Alb)  except  in  those  areas  influenced  by  local  topographic 
upwelling  (e.g.,  at  the  entrance  of  the  Bays)  or  in  particular  regions 
previously  vertically  mixed  by  the  passage  of  storms  (Montes-Hugo 
et  al.,  2005).  Likewise,  greater  phytoplankton  concentrations  as  a 
function  of  depth  were  commonly  observed  just  above  the  pycnocline 
depth  as  inferred  from  drastic  vertical  increase  of  CTD-derived  chi  o 
(>2  mgm  4 )  associated  with  major  depth  changes  on  sea  water  density 
(>0.01  kg  m  4)  (Appendix  A,  Fig  Alb).  Therefore,  these  environments 
offer  a  mosaic  of  case  studies  characterized  by  waters  with  different 
vertical  stratifications  and  consequently  different  vertical  structures  of 
optical  properties  and  biogeochemical  variables. 

22.  Aenol  surveys  ond  flight  mission  settings 

Airborne  Rn  (A,  0  +  )  and  5  (z)  measurements  were  obtained 
between  12:54  and  14:30  pm  local  time  and  over  waters  deeper 
(>50  m  depth)  than  penetration  depth  of  the  lidar  system  (Fig.  1).  To 
avoid  bottom-related  reflectance  and  lidar  backscattering  contribu¬ 
tions,  only  measurements  over  deep  waters  (i.e.,  >50  m  bottom 
depth)  were  analyzed.  Optimal  flight  weather  conditions  (i.e.,  cloud- 
free  skies,  wind  speed  <  4  m  $“  *)  were  checked  o  priori  based  on  local 
weather  forecast  (http://www.wunderground.com/)  to  maximize  the 
number  of  comparisons  between  passive  and  active  optical  measure¬ 
ments.  For  the  aerial  survey,  the  altitude  and  speed  were  300  m  and 
70  m  s  \  respectively.  Based  on  these  average  flight  characteristics, 
we  collected  a  total  of  3.0  x  104  passive  radiometric  measurements 
(upwelling  and  downwelling)  and  1.66  x  105  lidar  shots  along  a  total 
distance  of  332  km  during  83  min.  Airborne  optical  measurements 
were  geo-located  every  30  lidar  shots  (i.e.,  1  s)  during  the  full  survey. 

2.3.  Opticol  determinorions  oboord  the  otrp/one 

2.3.1.  Passive  measurements 

Upwelling  radiance  (Lu  (\))  and  downwelling  irradiance 
( Ed  (A)pUne)  were  measured  at  411,  443,  491,  509,  553,  665,  and 
780  nm  to  approximately  match  bands  1-7  of  the  spacebome  ocean 
color  sensor  SeaWiFS  (Churnside  &  Wilson.  2008).  The  sensors  (Micro 
Surface  Acquisition  System,  MicroSAS.  Satlantic)  have  bandwidths  of 
10  nm.  The  upwelling  radiance  sensor  has  a  typical  saturation  radiance 
ofSpWcm  2  nm  '  and  a  field  of  view  (FOV)  of  3°  (half  angle)  in  air. 
Based  on  this  FOV  and  a  sampling  rate  of  6  observations  per  second,  the 
pixel  diameter  of  Lu  MicroSAS  measurements  was  16  m  (along-track)  x 
200  m  (across  track).  The  MicroSAS  irradiance  sensor  has  a  typical 
saturation  irradiance  of  300  pW  cm  2  nm  1  and  a  noise  equivalent  of 
2.5x10  3pWcm  2nm  \  The  orientation  of  the  MicroSAS  sensors 
were  orthogonal  to  the  sea  surface  plane  (i.e.,  Lu  (X)  and  Fd  (A)^4™* 
were  measured  at  1 80°  and  0°  with  respect  to  zenith). 

2.32.  Active  measurements 

Lidar  backscattering  measurements  were  obtained  with  a  Fish 
Lidar  Oceanic  Experimental  (FLOE)  system  (Churnside  et  al.,  2001), 
looking  down  through  a  belly  port  of  a  twin-engine  aircraft.  FLOE 
detector  was  set  up  forward-looking  and  15°  off  vertical  to  minimize 
specular  reflections  from  the  sea  surface.  In  each  lidar  waveform,  the 
sea  surface  was  identified  based  on  the  large  increase  in  signal  when 
the  pulse  reached  the  surface.  In  this  study,  the  lidar  signature  at  the 
air-sea  interface  was  clearly  defined  due  to  the  absence  of  breaking- 
waves  and  fog  Hence,  after  eliminating  the  lidar  propagation  path 
in  air  (i.e.,  distance  between  airplane  and  sea  surface),  5  can  be 
computed  for  each  lidar  pulse  from  the  vertical  profiles  of  photo¬ 
cathode  current  by  applying  calibration  factors  related  to  optical 
system  parameters  (e.g.,  laser  pulse  energy,  surface  losses,  receiver 
area,  detector  responsivity).  The  FLOE  laser  is  linearly  polarized  and 
has  beam  divergence  and  receiver  field  of  view  of  1 7  mrad. 

The  receiver  is  polarized  in  the  orthogonal  plane  and  sampled  at  a 
rate  of  109  samples  per  second.  The  dynamic  range  of  the  receiver  is 
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sufficient  to  allow  5  measurements  to  a  depth  of  nearly  100  m  in  the 
clearest  ocean  waters.  In  our  coastal  study  area,  the  greater 
attenuation  of  laser  energy  reduced  FLOE  penetration  depth  to 
~30m  on  average.  The  green  (532  nm)  laser  source  was  pulsed 
with  an  energy  of  100  mj,  a  pulse  length  of  15  ns,  and  a  pulse 
repetition  frequency  of  30  Hz.  This  provides  a  pixel  size  of  5  m  with  an 
along-track  spacing  of  about  2  m.  The  sampling  vertical  resolution 
into  the  water  column  was  about  0.1  m. 

2.4.  Quality  control  of  oirbome  opticol  determ motions 

ln-situ  airborne  passive  optical  measurements  were  first  screened  to 
remove  data  affected  by  variable  illumination  caused  by  patchy  cloud 
layers  and  by  the  presence  of  whitecaps.  Additional  filtering  of  optical 
data  included  the  deletion  of  segments  characterized  by  airplane  turns 
(Fig.  1)  and  coinciding  with  drastic  spatial  changes  on  Ed  and 
magnitude  due  to  changes  on  orientation  of  airborne  sensors. 

Pixels  affected  by  clouds  were  screened  and  removed  from  further 
processing  by  comparing  field  Ed  (553)airp,ane  values  against  modeled 
Ed  (553)  above  the  sea  surface  (£d  (553)modd).  Ed  (553)modcl  estimates 
were  computed  based  on  a  radiative  transfer  model  (Hydrolight, 
Sequoia  Scientific,  Inc.)  using  local  meteorological  conditions  reported 
by  weather  stations  (visibility,  sea  level,  atmospheric  pressure,  wind 
speed,  relative  humidity,  Bob’s  P.C.  Connection,  Kodiak  Island,  Alaska, 
www.wunderground.com,  precipitable  water  content,  NOAA/GSD)  or 
extracted  from  satellite  imagery  (e.g.,  total  ozone  from  Total  Ozone 
Mapping  Spectrometer, TOMS)  during  the  initial  (21:54  GMT,  August 
17th)  and  final  (23:30  GMT,  August  17th)  stage  of  the  flight  survey. 
Solar  zenith  angle  was  calculated  using  a  solar  position  calculator  and 
based  on  date  and  location  (University  of  Oregon,  Solar  radiation 
monitoring  laboratory).  The  airmass  type  criterion  for  choosing 
aerosol  composition  was  based  on  wind  direction  data  provided  by 
Kodiak  Island  weather  station.  Aerosol  thickness  during  this  study 
was  very  small  (<0.04)  (SeaWiFS-derived  9  km  resolution,  wave¬ 
length  =865  nm)  and  is  consistent  with  the  lack  of  aerosol  contribu¬ 
tions  due  to  volcanic  plumes  during  this  study  (Alaska  Volcano 
Observatory,  http: //puff  images.aIaska.edu/puffAeries  airroutes/ 
Puff  20  years  1988  2008_series_conc.mpg). 

Glint  and  foam  contaminated  pixels  during  airborne  measure¬ 
ments  were  screened  based  on  Hydrolight  simulations  of  l ^  at  780  nm 
(Lu  (780)modd)  using  the  following  settings:  1 )  atmospheric  conditions 
similar  to  Ed  (553)modd  runs  but  wind  speed  equal  to  zero  (i.e.,  flat 
ocean),  and  2)  the  range  of  in  situ  chi  o  values  reported  in  this  area  (chi 
o  =  0.05  to  1 0  mg  m  3)  Magnitude  of  l ^  ( 780)  measured  outside  the 
theoretical  radiometric  interval  resulting  from  the  above  constraints 
was  considered  contaminated  with  glint  or  foam  contributions  and 
removed  from  further  processing.  Water  scattering  and  absorption 
coefficients  were  derived  from  Smith  and  Baker  (1981 ),  and  particulate 
scattering  and  absorption  coefficients  were  computed  from  the  Gordon 
and  Morel  empirical  equations  (Gordon  &  Morel,  1983.  Morel,  1991), 
and  interpolated  phytoplankton  specific  absorption  coefficients 
reported  by  Prieur  and  Sathyendranath  (1981).  Only  11.8%  of  the 
original  airborne  data  (i.e„  127  bins  with  250-m  resolution)  corre¬ 
sponded  to  the  ideal  approximation  ‘clear  skies  and  a  glint-foam  free 
ocean’  based  on  modeled  radiometric  thresholds  for  Ed  (553) 
(range  =  106.9  to  109.5 pW cm  2nm  ])  and  U,  (780)  (range  =  4.87 
10  2  to  8.86  10  2pWcm  2  nm  1  sr  !). 

2.5.  Atmospheric  correction 

A  quasi-single-scattering  approximation  (Rayleigh-aerosol  multi¬ 
ple  scattering  ignored)  was  used  to  relate  water-leaving  radiance  (Lw) 
to  cloud  and  glint-foam  free  iu  (Cordon  et  al.,  1983): 

4(M  =  4(M  +  4  +  t(\)Lw(\)  4-  L$int  (2) 


where  t(A)  is  the  diffuse  transmittance  of  the  atmosphere  below  the 
aircraft,  and  in  La  and  Lgh-„t  are  the  radiance  contributions  due  to 
Rayleigh/Fresnel,  aerosol,  and  glint,  respectively.  Lr  is  derived  from 
radiance  phase  functions  for  water  molecules  that  depend  on  incident 
solar  angles  (zenith  and  azimuth),  and  Fresnel  reflectance  estimates 
(Montes-Hugo  et  al.,  2005).  La  was  calculated  over  Clearwater  pixels 
where  a  minimum  water-leaving  radiance  at  665  nm  is  expected: 

La  =  W(\-,\0)(LU(665)-Lr(665))  (3) 

W(K|..\0)=£(\)F0(\)/F0(665)  (4) 

\o  =  665nm  and  A/ =  41 2,  443,  491,  509,  and  553  nm.  Assuming  a 
maritime  atmosphere,  £  is  1  (Gordon  et  al.,  1983).  F0(A)  is  the 
instantaneous  extraterrestrial  solar  irradiance  (ozone  optical  thick¬ 
ness  =  0)  (Neckel&  Labs,  1984). 

Glint  contribution  to  Lu  after  Hydrolight-based  pre-screening  of  Lu 
measurements  was  quantified  with  a  first-order  adjustment  by 
subtracting  Lu  (780)  from  (Lee  et  al.,  2001 ).  Skylight  path  radiance 
contribution  was  not  accurately  quantified  due  to  the  lack  of  sky¬ 
looking  radiance  sensor  but  it  was  assumed  small  due  to  the  relatively 
thin  atmospheric  layer  between  the  airplane  and  the  sea  surface.  In 
this  study,  ignoring  atmospheric  contributions  to  may  result  in  Lw 
overestimation  ofupto70%  in  shorter  wavelengths  (e.g.,  A  =  443  nm) 
due  mainly  to  molecular  backscattering,  and  48%  at  longer  wave¬ 
lengths  (e.g.,  A  =  553  nm)  because  of  the  larger  influence  of  aerosols. 

Assuming  a  negligible  attenuation  of  Ed  due  to  the  atmospheric 
path  below  the  airplane,  Rn  (A.  0  4 )  was  derived  as  L *  (A)  divided  by 
Ed  (A)atrplan*  (cosine  of  zenith  angle  =  1)  (Montes-Hugo  et  al.,  2005). 
In  case  1  waters,  Rn  (A,  0+ )  can  be  approximately  related  (~ 20%  bias) 
to  inherent  optical  properties  of  the  water  body  with  the  following 
expression  (Gordon  et  al .  1988): 

Ra(K,0J«0.54R/Q  (5) 

R/Q  =0.095  {hb{\)/(bb(\)  +  a(\))}  (6) 

where  bt  includes  light  backscattering  contributions  due  to  seawater 
and  particulates,  o  includes  light  absorption  contributions  due  to 
seawater,  colored  dissolved  organic  matter  and  particulates,  and  R/Q 
is  a  shape  distribution  factor  that  is  influenced  by  the  light  field 
geometry. 

2.6.  Comparisons  of  optica/  parameters  derived  from  passive  o  rid  active 
airhame  sensors 

One  important  goal  of  this  study  was  to  understand  how 
information  provided  by  passive  and  active  optical  sensors  was 
related  despite  the  differences  in  spatial  (i.e.,  vertically  integrated 
from  MicroSAS  versus  vertically-resolved  from  FLOE)  and  spectral 
(e.g.,  multi-wavelength  MicroSAS  versus  single-wavelength  FLOE) 
resolution. 

Remote  sensing  reflectance  is  an  apparent  optical  property  (AOP) 
that  is  affected  by  light  field  angular  characteristics  such  as  solar 
zenith  angle.  Therefore,  a  more  direct  way  to  evaluate  the  effects  of 
water  characteristics  on  MicroSAS-denved  Rn  and  FLOE-derived  S 
measurements  independently  from  illumination  effects  is  to  use  Rn 
spectral  bands  and  curvature  ratios  (Montes-Hugo  et  al.,  2005).  A 
preliminary  sensitivity  analysis  suggested  two  main  MicroSAS- 
derived  variables:  ^(443,  0+)/R„(490,  04)  (hereafter  R\)  and  Rn 
(508,  0-1- )//?„(  553,  0+)  (hereafter  R2)-  Interestingly.  Rt  is  inversely 
related  to  the  average  size  of  particulates  of  a  specific  water  parcel 
(Gordon  &  Morel,  1983)  while  R2  is  inversely  related  to  the 
concentration  of  pigmented  particulates  containing  chlorophyll  o 
(Montes-Hugo  et  al.,  2005). 

The  comparisons  between  passive  and  active  optical  measure¬ 
ments  were  organized  according  to  three  main  questions:  How  do  the 
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magnitude  of  MicroSAS  and  FLOE  optical  parameters  vary?  (1),  How 
does  a  relate  to  MicroSAS- derived  lOPs?  (II),  What  is  the  depth  of  the 
lidar  scattering  layer  having  the  largest  influence  on  the  magnitude 
and  spatial  variability  of  Micro SAS-de rived  lOPs  and  Rn  ratios?  (111). 

Regarding  (1),  we  calculated  the  moments  around  the  mean 
(arithmetic  average,  standard  deviation,  skewness  and  kurtosis)  of 
the  Rn  ratios  and  S.  Skewness  (t/0  and  kurtosis  (77)  represent  the  third 
and  fourth  standardized  moments  (m),  respectively: 

4 >  =  ^/(J3  (7) 

il  =  m/<i4  (8) 

Mrn  =  J  (x-n)m/(x)dx  (9) 

where  is  the  arithmetic  average,  a  is  the  standard  deviation,  x  is  the 
random  variable,  and  f[x)  is  the  continuous  univariate  probability 
density  function.  Large  1/7  relative  to  a  Gaussian  curve  (1/7  =  0) 
corresponds  with  an  ’excess’  of  signal  that  may  be  related  to  a  greater 
concentration  or  composition  of  optical  components  (e.gM  fish  species 
and  abundance  in  lidar  waveforms).  Likewise,  high  77  with  respect  to  a 
Gaussian  function  (77  =  3)  means  that  a  greater  proportion  of  the 
variance  is  concentrated  around  the  mean,  and  implies  a  more 
’homogeneous’  origin  of  the  optical  signal. 

For  each  lidar  profile,  a  was  derived  from  the  regression  slope  of 
neperian  log-tra  ns  formed  S  as  function  of  depth  after  automatic  and 
manual  deletion  of  those  1  -m  vertically  integrated  lidar  returns  where 
S  increases  with  depth.  In  case  where  S  increased  with  depth  a  was 
not  calculated.  Automatic  pre-screening  and  QC  of  S  waveforms  was 
made  by  setting  flags  every  1-m  depth  and  detecting  positive 
differences  of  S  with  depth  as  the  lidar  signal  propagates  deeper  in 
the  water  column. 

The  horizontal  patchiness  of  S  (hereafter  k)  was  derived  from  the 
slope  of  the  log-transform  power  spectrum  of  vertically  averaged  S 
(from  2  to  20  m  depth)  as  a  function  the  log-tra  ns  formed  spatial 
frequency.  The  coefficient  of  variation  of  S  (cv)  was  computed  as  the 
standard  deviation  of  the  signal  without  spatial  integration  divided  by 
the  arithmetic  average  ofS.S  parameters  (i.e.,  kurtosis,  skewness,  «,  k, 
vertically  integrated  S,  and  cv)  were  computed  using  the  full  FLOE 
vertical  resolution  (i.e.,  179  data  points  per  profile)  and  between  2 
and  20  m  depth  to  minimize  surface  effects  caused  by  air  bubbles  and 
large  noise  contribution  beyond  lidar  penetration  depth  (i.e.,  ~1Q 
standard  deviation  above  the  lidar  noise). 

To  examine  (11),  bb  and  a  coefficients  were  calculated  from 
MicroSAS-derived  Rn  at  the  lidar  wavelength.  Only  MicroSAS  data 
corresponding  to  lidar  profiles  with  monotonic  decrease  of  S  with 
depth  were  selected.  An  updated  version  of  the  quasi-analytical 
algorithm  (QAA  v5)  was  used  to  convert  Rn  into  lOPs  (Lee  et  al.. 
2009).  Briefly,  QAA,v5  protocol  has  the  following  calculation 
modules:  1)  Rn  just  below  the  sea  surface,  2)  spectral  shape  of 
and  particulate  backscattering  (6bp  (\)),  3)  spectral  shape  of  a  (\),and 
4)  partition  of  a  (\)  in  detritus  and  phytoplankton  contributions. 
Based  on  the  NOMAD  validation  dataset,  QAA  v5  can  retrieve  a’s  in 
the  visible  range  with  a  coefficient  of  determination  greater  than  0.9 
and  root-mean  square  difference  (log  scale)  smaller  than  0.2  (Lee 
et  al ,  2009).  Once  bb  (532)  and  a  (532)  are  derived,  Kd  (532)  can  be 
easily  estimated  for  cloud-free  skies  (Gordon,  1989). 

To  answer  (111),  we  performed  two  different  analyses.  First, 
MicroSAS-derived  lOPs  were  related  to  smoothed  S  curves  using 
vertically  integrated  1-m  bins  (52$  (zmax)  =  sum  of  the  product 
Sxdz  from  vertical  layer  i  to  i  +  DZ,  where  zmax.  dz  and  DZ  are  the 
upper  limit  of  the  integral,  the  depth  differential,  and  the  smoothing 
interval  =  1  m,  respectively).  Lastly,  the  influence  of  different  vertical 
shapes  ofS  (z)  on  MicroSAS-derived  Rn  was  examined  based  on  the 


following  contrasting  cases:  1)  ’exponential’  versus  ’non-exponential’ 
decrease  of  51 S  (zmax)  with  depth,  2)  ’deep’  (i.e.,  ^  11  m)  versus 
’shallow’  52$  (zmax)  maximum,  3)  ’thick’  (i.e.,^  5  m)  versus  ’thin’  52$ 
(zmax)  maximum,  and  4)  ’strong’  (i.e.,  ^7  10  4  m  sr  ')  versus  ’weak’ 
51 S  ( zmax )  maximum.  These  binary  categories  are  defined  with  respect 
to  the  median  of  the  sample  and  considering  profiles  with  only  one  52  S 
(zmax)  maximum  peak  between  3  and  20  m  in  depth.  The  magnitude  of 
each  5 ZS  (zmax)  maximum  (i.e,  largest  peak  height)  is  relative  and 
defined  with  respect  to  a  baseline  delineated  by  the  exponential  decay 
of  the  background  signal.  The  width  of  each  52!  5  (zmax)  maximum  is 
estimated  at  the  base  of  the  Gaussian-type  peak  and  corresponds  to  the 
distance  between  two  intersecting  points  made  by  the  peak  tails  and  the 
background  signal  Notice  that  the  final  number  of  lidar  profiles  with 
different  types  of  55  S  (zmax)  maximum  was  relatively  small  due  to  the 
reduced  number  of  coincident  FLOE-MicroSAS  measurements  without 
atmospheric  interference  (-10%  of  original  dataset).  Although  FLOE 
waveforms  with  a  dominant  negative  exponential  term  were  observed 
in  numerous  occasions  (e.g,  >500  bins),  it  was  decided  to  use  a  number 
of  observations  comparable  to  those  obtained  in  the  ’non-exponential’ 
5 ZS  (zmax)  case  study  to  achieve  an  even  weighting  during  statistical 
comparisons.  Post-comparisons  of  moments  (low  and  high)  around  the 
mean  between  these  lidar  profile  types  and  using  all  ’exponential’ 
samples  did  not  change  observed  differences  between  exponential  and 
non-exponential  lidar  waveforms 

2.7.  Statistical  analysis 

Passive  and  active  airborne  spectral  optical  products  were 
aggregated  into  250-m  horizontal  bins  (i.e,  125  FLOE  profiles,  22 
MicroSAS  measurements)  before  performing  spatial  correlations  and 
variance  decomposition.  The  choice  of  250-m  bins  was  made  for 
several  reasons:  1)  It  is  the  best  resolution  of  satellite  ocean  color 
sensors  with  relatively  short  (1  to  several  days)  revisiting  time  (e.g, 
MOD1S),  2)  The  sample  size  per  bin  is  relatively  large  (FLOE: 
N=  125x200,  MicroSAS:  N=360)  for  histograms  and  statistical 
tests,  and  3)  Smaller  bin  sizes  would  be  more  influenced  by  small- 
scale  surface  effects  (e.g,  glint,  bubbles  due  to  wave  breaking), 
producing  more  noise/signal  since  these  surface  contributions  cannot 
be  completely  corrected.  The  sensitivity  of  the  along-track  variability 
of  lOPs  to  horizontal  changes  of  vertically  integrated  S  (i.e,  525 
(zmax))  was  quantified  using  principal  component  analysis  (PCA) 
(Pearson,  1901).  The  contribution  of  each  variable  to  the  first  3 
orthogonal  axes  (i.e,  those  who  explain  >80%  of  total  variance)  was 
interpreted  based  on  the  loading  of  each  variable  to  each  eigenvector 
(i.e,  principal  component  or  PC).  Fraction  of  total  variability  explained 
by  principal  components  decreases  from  PCI  to  PC3.  For  each 
eigenvector,  the  sign  of  the  loading  of  each  optical  parameter  was 
used  to  establish  the  sign  of  the  correlation  (i.e,  direct  or  inverse) 
between  MicroSAS  and  FLOE  variables.  The  intensity  of  the  relation¬ 
ships  between  MicroSAS  and  FLOE  optical  variables  were  measured 
using  Spearman  Rank  non- para  metric  correlation  coefficients  (ps) 
(Spearman.  1904).  Comparisons  were  only  performed  using  a  spatial 
lag  equal  to  zero  due  to  lidar  data  gaps  (e.g,  FLOE  profiles  between 
21:20  and  22:33  GMT)  along  the  flight  track.  Correlations  were  also 
examined  between  different  characteristics  (amplitude,  g,t  depth,  g2, 
and  width,  g3)  of  S  waveforms  having  a  single  sub-surface  peak. 
Assuming  a  Gaussian  mode!  (C).  the  relationship  between  gl,  g2  and 
g3  is: 

C(z)  =  gl(?-05«*-*'/*>\  (10) 


The  amplitude  coefficient  is  sometimes  expressed  as  h/(2u  g3), 
where  h  is  the  total  vertically  integrated  value  of  the  property  above  a 
baseline  (Millan-Nufiez  et  al,  1997). 
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The  significance  of  change  in  the  averaged  spectral  ratios  and  lidar- 
derived  parameters  with  different  vertical  S  distributions  (see 
Section  2.6)  was  evaluated  at  a  confidence  interval  of  probability  of 
95%  and  99%  and  based  on  t-Student  comparisons  with  1  -tail.  The  null 
hypothesis  was  M<groUp  i)  =  M<group  2).  with  group  l,2...n  equal  to 
different  subsets  of  optical  data  corresponding  to  pre-defined  S  (z) 
clusters.  Change  of  skewness  and  kurtosis  in  each  dataset  was 
measured  with  respect  to  a  Gaussian  distribution,  i/r  or  77  were 
significant  with  a  5%  error  (Type  1)  when  their  absolute  magnitude 
was  greater  than  two  standard  errors  of  t/r  (ses  =  y/6 /  N,  where 
N  =  number  of  observations)  and  T]  (sek  =  y/24 /  N),  respectively 
(Tabachmck  &  Fidell,  1996).  Differences  in  the  coefficient  of  variation 
of  S  per  bin  between  two  case  studies  with  N  observations  ( i.e.,  N 1  and 
N2)  were  examined  by  comparing  a  pair  of  proportions  (i.e.,  pi  and 
p2)  with  a  one-side  t-test  (|t|  =  V^N1N2  /  (N 1  +  N2)]pl  — p2|  /  ^fq), 
p  =  (p]  N 1  +p2  N2)/(N1  +  N2),  q  —  1  —p,  degrees  of  freedom  =  N1  + 
N2-2  (STAT1STICA  software). 


scattering  layers  were  discontinuous  (e.g.,  Lidar  returns  at  1.23  and 
3.1  km  with  respect  to  the  starting  distance  of  the  lidargram)  (Fig  2a, 
d).  Conversely,  lower  R\  values  consistently  followed  deepening  of  a 
mid-water  (5  to  15  m  depth)  scattering  layer  with  an  S  magnitude 
comparable  to  surface  values  (i.e.,  >1  10  5  m  1  sr  ]).  At  a  coarser 
horizontal  spatial  resolution  (i.e.,  250  m),  the  median  and  standard 
deviation  of  Rt  still  evidenced  the  presence  of  this  submerged 
scattering  feature  (Fig.  2b).  This  high  reflective  ’thin’ 
layer  (thickness ~  1.5  m)  became  shallower  at  the  end  of  the  lidar 
sampling  segment  (i.e.,  3.7  to  4.1  km)  and  caused  a  drastic  increase 
of  Rt  sd  (up  to  33-fold)  and  cv  (up  to  24-fold).  In  contrast,  minimum 
R,  sd  (<5  10  3)  and  cv  (<0.5%)  were  associated  with  lidar  profiles 
having  a  monotonic  decrease  of  S  with  depth  (0  to  0.5  km)  or 
relatively  weak  (~  1.5  10  7  m  1  sr  ')  and  deep  (14-16  m  depth) 
scattering  layers  (e.g.,  2  to  2.5  km)  Interestingly,  maximum  values  of 
higher  moments  around  R,  mean  evidenced  a  horizontal  shift  of  250 


3.  Results 

3.1  Correspondence  between  magnitude  ond  spotiol  vonobility  of 
spec  fro/  reflectance  rotios  ond  lidor  bockscottering  signol 

Based  on  spectral  ratios  of  remote  sensing  reflectance,  the  color  of 
the  ocean  was  correlated  with  the  magnitude  of  lidar  signal  extinction 
and  degree  of  S  patchiness  along  the  horizontal  component  (Table  1 ). 
A  typical  example  showing  the  coupling  between  passive  and  active 
optical  information  at  different  horizontal  spatial  scales  ( 16  to  250  m) 
is  illustrated  in  Fig.  2.  The  airborne  measurements  were  obtained 
nearby  Whale  Island  (58.00°  N,  -  152.74°  W,  Fig.  1)  at  22:36  h  GMT 
and  over  relatively  deep  waters  (bottom  depth  =  120  m).  In  general, 
R,  had  a  closer  association  with  fine  structure  changes  on  2-D 
distributions  ofS  (i.e.,  lidargrams)  than  did  R2  (Fig  2a).  Likewise,  the 
highest  Rt  values  commonly  coincided  with  lidar  shots  characterized 
by  drastic  attenuation  with  depth  and  where  surface  (0-5  m  depth) 


Table  1 

Correlation  (/>s)  between  airborne  passive  and  active  optical  parameters.  Definitions  of 
I,  II  and  III  comparisons  are  defined  in  Section  2.6  of  Methods.  pt  significant  at  95% 
(P<  0.05)  and  99%  (P< 0.01 )  confidence  level. 
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Fig.  Z  Spatial  relationship  between  airborne  remote  sensing  reflectance  ratios 
(R,  and  R2)  and  vertical  distribution  of  lidar  back  scattering  along  the  flight  track, 
a)  Rn  ratios  calculated  every  16  m,  b)  ft  and  sd  of  Rn  ratios  calculated  every  250  m. 
c)  same  as  b)  but  for  kurtosis  and  skewness.  Gaussian  kurtosis  (dotted  line)  and 
skewness  (solid  line),  d)  lidar-derived  S  profiles,  S  in  neperian  log  scale,  distance 
between  profiles  is  2  m. 
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Table  2 

Principal  component  analysis  between  passive  and  active  optical  measurements.  In 
bold  are  indicated  those  variables  with  the  largest  contribution  to  each  eigenvector 
(i.e.,  *|0.3|. 


Sensor 

Variables 

Histogram  metrics 

PCI 

pa 

pa 

MicroSAS 

ft* (532) 

X 

0.16 

039 

-035 

a (532) 

X 

0.11 

0.41 

-035 

FLOE 

ES(3) 

0.26 

-033 

-0.26 

ES(5) 

0.20 

-035 

-037 

ES(I0) 

-0.46 

-0.08 

-030 

ES(15) 

-0.44 

-0.16 

-038 

EStfO) 

-038 

-0.21 

-0.14 

to  500  m  with  respect  to  median  and  sd  (Fig.  2c).  Roughly,  the 
largest  kurtosis  and  skewness  of  tended  to  lead  peaks  on  Rj  median 
and  sd  but  this  was  not  always  the  case  (e.g.,  r)  =  4.6,  tf/  =  1 .42,  median 
R,  =  1.06,  std  Ri  =4  10  \  first  bin  centered  at  0.41  km). 

R,-a  relationships  suggested  a  strong  association  between  waters 
dominated  by  small-sized  particulates  and  greater  vertical  attenua¬ 
tion  of  the  lidar  signal  (p3  =  0.38).  As  inferred  from  k,  greater 
horizontal  patchiness  of  lidar  backscattering  was  associated  with 
water  containing  larger  particles  ( ps  for  R,  —  k=  -0.31 ).  The  total 
strength  of  the  lidar  signal  integrated  from  2.1  to  20  m  was  not 
significantly  correlated  with  any  spectral  Rn  ratio.  However,  the 
magnitude  of  appeared  to  be  related  to  the  horizontal 
homogeneity  of  R!  and  R2  within  each  bin  (Spearman  rank  correla¬ 
tions  between  standard  deviation  of  Rn  ratios  and  £  S  up  to  0.36).  No 
spatial  coherence  was  detected  between  magnitude  of  MicroSAS- 
derived  lOPs  and  lidar  attenuation  coefficient  values  (P>0.05, 
Table  1).  Neither  were  magnitudes  of  R„-based  lOPs  explained  by 
the  magnitude  of  differences  on  S  at  depths  above  20  m.  Inverse 
relationships  between  MicroSAS  and  FLOE  optical  signatures  were 
clearly  manifested  at  20  m  depth  where  correlations  between  bb 
(532),  a  (532)  and  S  were  negative  (comparison  Hi  in  Table  1 ). 

In  this  study  we  applied  PCA  to  identify  major  sources  of  spatial 
variability  affecting  concurrent  FLOE  and  MicroSAS  measurements 
(Table  2 ).  The  first  Principal  component  (i.e.,  PCI )  was  mainly  affected 
by  variability  of  the  'deep*  (10-20  m  depth)  lidar  backscattering. 
Unlike  PCI,  PC2  and  PC3  were  mainly  affected  by  variations  on 
radiometer-derived  lOPs  and  vertically  integrated  lidar  volume 
backscattering  near  the  surface  (3  to  5  m).  In  PC2,  an  inverse 
relationship  between  passive  ( bb  (532),  a  (532))  and  active  (£S  at  3 
and  5  m  depth)  parameters  suggests  that  horizontal  changes  of  lOPs 
were  not  associated  with  inter-bin  changes  of  lidar  attenuation  (i.e., 
a).  PC3  revealed  the  linkage  between  of  MicroSAS-derived  lOPs  and 
lidar  volume  backscattering  at  intermediate  depths  (^  3  and  <  10  m). 


Interestingly,  bb  (532)  and  a  (532)  had  comparable  contributions  to 
PC3,  which  indicates  a  modulation  of  lidar  signal  variability  due  to 
spatial  fluctuations  of  both  backscattering  and  absorption  coefficients. 

32.  Changes  on  ocean  co/or  statistical  moments  due  to  dijfferenr  ’shapes’ 
on  lidar  waveforms 

Spatial  patterns  of  R„  ratios  measured  in  this  study  differ 
depending  on  the  number  of  lidar  derived  backscattering  layers 
with  depth.  One  single  layer  with  a  homogeneous  composition  of 
optical  constituents  (i.e.,  S  only  dominated  by  one  exponential  decay 
term)  extending  from  the  sea  surface  down  to  20  m  typically 
corresponded  to  averaged  higher  R,  (t=  5.95,  P<0.01 ),  and  higher 

kurtosis  and  skewness  of  R,  (tj  =  3.20.  \p  =  0.19)  than  a  vertical  multi¬ 
layered  system  of  S  where  the  amplitude  of  lidar  signal  varied  with 
depth  (Table  3,  Fig.  2c-d).  Waters  without  sub-surface  S  'bumps'  also 
had  on  average  2-fold  lower  integrated  S  over  the  lidar  penetration 
range  (t  =  5.3,  P<0.01).  In  general,  lidar  profiles  with  more  complex 
vertical  structure  produced  S  with  larger  coefficient  of  variability  (cv) 
(t  =  8.95,  P<0.01)  and  higher  k  values  (t=  10.46,  P<0.01 )  (Table  3). 

In  general,  locations  with  deeper  S  peaks  (smoothed  to  l~m 
vertical  x  250  m  horizontal  resolution)  had  lower  horizontal  disper¬ 
sion  (i.e.,  cv)  (t  =  3.1 6,  P<0,01)  of  targets  contributing  to  S  (Table  4). 
These  S  distributions  tended  to  be  more  leptokurtic  and  symmetric, 
and  commonly  associated  with  larger  R,  kurtosis  and  negative  R, 
skewness.  The  group  of  lidar  profiles  classified  by  'weak*  sub-surface 
maxima  shared  various  characteristics  of  the  'deep'  cluster.  Within  the 
’weak’  lidar  peak  category,  R,  histograms  had  lower  skewness  (t  = 

—  2.58,  P  =  0.02),  the  total  integrated  S  between  2.1  and  3  m  was  62% 
smaller  (t=  — 3.0,  P<0.01),  and  the  5  distribution  was  less  patchy 
(t=  - 1.85,  P=0.04)  and  more  homogenous  (i.e.,  lower  cv)  horizon¬ 
tally  (t=  —4.90,  P=0.03)  than  the  'strong'  cluster.  For  those  lidar 
profiles  with  a  sub-surface  S  peak  relatively  narrow  ('thin'  scenario), 
Ri  was  higher  (t  —  2.0,  P=0.03),  R2  was  lower  (t=  1.8,  P=0.045),  S 
vertically  integrated  per  bin  tended  to  be  lower  (-30%).  S  skewness 
was  close  to  zero,  and  horizontal  distribution  of  S  was  less  patchy  than 
the  'thick'  type  lidar  peak  cluster.  Overall,  depth  and  magnitude  of  S 
sub-surface  peaks  were  inversely  correlated  (ps=  —  0.48,  t(N  — 2)  = 

—  2.52,  P=  0.02)  (Appendix  A.  Table  Al ). 

4.  Discussion 

Finding  relationships  between  optical  properties  derived  from  R„ 
(A,  04- )  and  lidar  profiles  is  an  important  step  forward  to  elucidate  1 ) 
the  effects  of  vertical  structure  on  ocean  color  imagery  and  satellite 
derived  products,  and  2)  the  composition  of  the  lidar  signal.  The  main 
results  of  this  contribution  were  organized  in  two  sections:  A)  How 


Tabic  3 

Airborne  reflectance  derived  parameters  over  waters  with  A(\  relatively  constant  (HS)  or  varying  with  depth  (IS):  mm  and  max  are  minimum  and  maximum  values,  respectively. 
R)t  R2.  k,  and  cv  (dimensionless).  bb{532).  a(532)  and  Krf(532)  (m  1 ).  (m  sr  1 )  10  ,.o(m'1). 


Ki 

r2 

(532) 

a  (532) 

Ka  (532) 

ES 

cr 

K 

cv 

HS 

min 

1.01 

0.68 

0.003 

0,06 

0.09 

352 

0.13 

0.03 

0.03 

max 

139 

233 

0.011 

0.60 

037 

7.08 

023 

137 

0.18 

n 

41 

41 

41 

41 

41 

41 

41 

it 

1,18 

1.16 

0.006 

0.16 

024 

5.08 

0.18 

0.56 

0.07 

T? 

3.20  ±1.87 

2.76  ±1.46 

7.74  ±4.51 

* 

0.19  ±036 

0.10±  127 

023  ±0.59 

IS 

min 

0.96 

1.07 

416 

0.62 

0.05 

max 

1.14 

133 

24.17 

221 

1.62 

n 

23 

23 

23 

23 

23 

23 

23 

23 

* 

1.04 

122 

10.11 

1.51 

0.58 

tl 

037  ±0.68 

2.43  ±033 

6  59  ±6.48 

* 

—  031  ±131 

-0.02  ±0  45 

0.64  ±0.96 

2S90 
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Table  4 

Airborne  reflectance  derived  parameters  over  waters  with  typical  'shallow',  'deep  .  “weak’. 'strong*. 'thin'  and  'thick'  lidar  profile  categories  (Section  2.6)  Between  parentheses  is  the 
number  of  bins  used  in  each  lidar  profile  category. 


Category 

Histogram  metrics 

MicroSAS  variables 

Histogram  metrics 

FLOE  variables 

•Shallow'd!) 

X±sd 

Ri 

1.05  ±0.06 

*±sd 

Es 

10.8  ±6.4 

fj±sd 

*1 

0.1S±0.81 

fj±sd 

s 

3.03  ±4. 11 

*±sd 

R| 

031  ±0.88 

■i±sd 

s 

0.82±  1.15 

*±sd 

R2 

130  ±0.05 

X±sd 

K 

7.05  ±11.00 

rj±sd 

r2 

238  ±0.62 

X±sd 

cv 

039 ±033 

*±sd 

r2 

-0.07  ±043 

’Deep’ (12) 

X±sd 

*1 

1.04  ±0.05 

X±sd 

ES 

93  ±  13 

rj±sd 

Ri 

037  ±03S 

rj±sd 

s 

9.87  ±  6.6S 

</r±sd 

Ri 

-0.69  ±131 

*±sd 

s 

0.48  ±0.78 

X±sd 

R2 

133  ±0.07 

X±sd 

K 

130  ±0.44 

rj±sd 

R2 

2.48  ±1.00 

X±sd 

cv 

0.37  ±032 

*±sd 

r2 

0.03  ±0.49 

Weak’  (10) 

X±sd 

Ri 

1.05  ±0.06 

Jt±sd 

ES 

73  ±1.7 

rj±sd 

Ri 

038  ±0.64 

r]±sd 

s 

8.69  ±733 

*±sd 

-0.87  ±1.17 

*±sd 

s 

030  ±0.81 

X±sd 

*2 

131  ±0.08 

X±sd 

K 

130  ±0.49 

fj±sd 

2.36  ±1.11 

X±sd 

cv 

039 ±036 

*±sd 

R 2 

-0.05  ±0.44 

X±sd 

•Strong’  (13) 

X±sd 

R. 

104  ±0.05 

X±sd 

Es 

12.1  ±5.1 

rj±sd 

R< 

036  ±0.73 

f]±sd 

s 

4  98  ±  S30 

*±sd 

Rj 

030  ±1.00 

*±sd 

s 

0.99  ±0.96 

X±sd 

R  2 

132  ±0.05 

X±sd 

K 

4.98  ±530 

rj±sd 

R  2 

2.49  ±0.56 

X±sd 

cv 

0.81  ±0.41 

*±sd 

R  2 

0.01  ±0.49 

X±sd 

Thin’  (9) 

X±sd 

R> 

1.08  ±0.06 

X±sd 

Es 

6.8  ±4.6 

rj±sd 

R. 

0.01  ±0.71 

r?±sd 

s 

S  03  ±  S.04 

*±sd 

R. 

-0.49  ±137 

i^±sd 

s 

0.10  ±1.09 

*±sd 

R  2 

1.18±0.06 

X±sd 

K 

136  ±6.49 

Tf±Sd 

R  2 

234  ±0.64 

X±sd 

cv 

0.70  ±034 

*±sd 

R  2 

-0.07  ±036 

Thick-  (14) 

*±sd 

R, 

1.03  ±0.04 

X±sd 

Es 

11.1  ±4.S 

rj±sd 

R. 

0.44  ±0.63 

r|±sd 

s 

839  ±535 

sd 

R. 

—  0.17±  1.02 

tjf±s  d 

s 

0.71  ±0.94 

X±sd 

R  2 

133  ±0.06 

X±sd 

K 

4.41  ±9.31 

rj±sd 

R2 

233  ±0.96 

X±$d 

cv 

030  ±0.28 

0±$d 

R2 

0.01  ±  0.S2 

strong  is  the  relationship  between  total  lidar  volume  backscattering 
(magnitude  and  spatial  variability)  and  lOPs  (e.g.,  bb  and  a 
coefficients)  as  estimated  by  passive  optical  systems?,  and  B)  Can 
we  obtain  information  about  vertical  shape  of  profiles  based  on  higher 
statistical  moments  of  airborne-derived  Rn  (A,  0  +  )  (e.g.,  skewness, 
kurtosis)? 

4.1 .  Spafio/  changes  on  lOPs  and  vertical  attenuation  of  lidar 
backscattering 

The  two-way  propagation  of  lidar  signal  through  a  finite  aqueous 
optical  path  involves  attenuation  of  photons  due  to  particulate  and 
dissolved  components  (see  exponential  term  in  Eq.  (1)).  Our  results 
suggested  that  MicroSAS-derived  lOPs  (a  (532)  and  bb  (532) 
coefficients,  and  Kd  (532))  were  weakly  correlated  ( ps  up  to  T0.15, 
p=  0.10)  with  lidar  attenuation  coefficient.  Based  on  the  FLOE  field  of 
view,  a  can  be  approximated  by  Kd  (532)  (Chumside  et  al.,  1 998),  thus 
a  is  expected  to  be  primarily  influenced  by  absorption  of  the  lidar 
signal  with  depth  due  to  background  optical  constituents  (e.g., 
phytoplankton).  That  was  not  the  case  as  R2,  an  index  of  phytoplank¬ 
ton  pigment  concentration,  was  uncorrelated  with  a  (Table  2). 
Therefore,  a  greater  phytoplankton  absorption  of  green  light  (i.e., 
wavelength  =  532  nm)  did  not  appear  to  be  a  major  factor  influencing 
the  lidar  signal  extinction.  However,  based  on  principal  component 
analysis  (PC3,  Table  2),  it  can  be  said  that  spatial  variability  of  a,  bb, 
and  Kd  derived  at  a  single  wavelength  covaned  with  each  other  but 
there  were  other  factors  affecting  these  relationships  in  the  first  two 
orthogonal  PCA  axis.  These  additional  factors  may  include  differences 
in  polarization  between  sensors  (e.g.,  cross  polarized  in  FLOE,  and 


unpolarized  in  MicroSAS)  and  phase  function  variations  caused  by 
different  illumination  geometries  (e.g.,  lidar  beam  normal  to  the  sea 
surface  versus  zenith  angle  of  sunlight  rays).  Also,  some  assumptions 
like  vertical  homogeneity  of  10P  distributions  and  a  constant  A(i  (z) 
may  not  be  valid  in  some  cases,  and  may  introduce  bias  on  Kd  and  a 
estimates  due  to  different  vertical  structure  of  optical  properties. 
More  research  will  clarify  based  on  our  optical  system  whether  the 
observed  poor  Kd-a  dependency  is  related  to  the  presence  of  lOPs 
inhomogeneities  and  A(i  changes  along  the  vertical  component  and/or 
to  variations  inherent  to  each  optical  sensor. 

Despite  the  relatively  weak  connection  between  magnitude  of  Rrs- 
denved  lOPs  and  integrated  lidar  backscattering,  higher  a  or  £ S 
values  were  associated  with  a  large  variability  of  R]  and  R2  within  a 
bin  (Table  1 )  In  NGOA  waters  and  during  the  same  period  of  the  year, 
Brown  et  al  (2002)  concluded  that  spatial  variability  of  lidar  signal 
appeared  to  increase  with  zooplankton  density  due  to  a  higher  degree 
of  spatial  patchiness  and  variation  in  zooplankton  density  within 
patches.  Greater  spatial  variability  on  zooplankton  distributions  has 
been  reported  in  regions  characterized  by  patchy  distribution  of  food 
(e.g.,  phytoplankton)  (Folt  &  Burns,  1999).  Therefore,  a  larger  spatial 
heterogeneity  of  phytoplankton  and  other  particulates  affecting  Rn 
ratios  can  be  expected  in  those  waters  characterized  by  more 
heterogeneous  spatial  distribution  patterns  of  zooplankton.  The 
impact  of  zooplankton  (i.e.,  an  important  optical  component 
affecting  S)  on  phytoplankton  (i.e.,  an  important  optical  component 
affecting  Rn  ratios)  was  presumably  evidenced  by  principal  compo¬ 
nent  analysis.  The  inverse  relationship  between  lOPs  and  ]TS 
eigenvector’s  loadings  in  PC2  was  probably  indicative  of  phytoplank¬ 
ton  depletion  and  consequently  a  reduction  due  to  zooplankton 


MA  Montes- Hugo  et  ol.  /  Remote  Sensing  of  Environment  J 14  (2010)  2584-2593 


259] 


0  001  0.0  !  0.001  0.01 


y.S'Unux)  (m  sr*1) 


]T.V(zituixi  (in  sr~! ) 


Fig.  3.  Classification  of  lidar  waveforms  with  a  single  sub-surface  S  peak,  a)  ‘Exponential’  versus  ‘non-exponential*  decrease  of  S  with  depth,  b)  shallow’  versus  ’deep',  c)  'weak'  versus 
‘strong’,  and  d)  ’thin’  versus  ‘thick*,  ji  (solid  line)  and  ±2  standard  errors  (se)  (dotted  line)  of  each  S  profile  smoothed  every  I  m.  Between  parentheses  is  S  ±  se  of  each  peak  Thex- 
axis  is  in  log  scale  with  base  10  and  £S  (zmax)  values  were  offset  in  +0.001  to  plot  *  —  2  se<0. 


grazing.  £S  calculated  at  3  and  5  m  depths  suggested  that  this  effect 
was  in  part  propelled  by  feeding  activity  of  surface  zooplankton  layer 
between  dawn  and  dusk.  However  during  the  day,  vertical  migrating 
zooplankton  (e.g.,  copepods)  like  those  inhibiting  NCOA  are  prefer¬ 
entially  situated  in  deeper  waters  (negative  loadings  of  ]£S  to  PCI 
between  10  and  20  m  in  depth).  Thus,  a  major  factor  defining  lOPs 
decrease  due  to  XI 5  increase  was  likely  the  feeding  activity  of 
zooplankton  during  the  previous  night  when  a  larger  volume  of 
grazers  was  closer  to  the  ocean  surface.  Spatial  variations  of  the 
production  of  pigmented  particulates  and  dissolved  colored  com¬ 
pounds  affecting  a  and  b*,  coefficients  did  not  appear  to  be  a  major 
mechanisms  explaining  magnitude  of  a  or  £ 5  values  since  direct 
relationships  between  these  parameters  were  secondary  as  evidenced 
in  PC3. 

4.2.  Response  of  spectra  reflectance  ratios  to  variation  on  vertical  'shape' 
of  Ildar  waveforms 

Spatial  patterns  of  ocean  color  in  our  study  area  changed  if  sub¬ 
surface  exponential  decline  ofS  as  a  function  of  depth  was  interrupted 
by  high  scattering  layers.  This  phenomenon  was  particularly  intrigu¬ 
ing  when  analyzing  the  horizontal  shift  between  maximum  of  low 


(median  and  sd)  and  high  (?/and  i//)  moments  of  R j  (fig.  2c-d)  It  is 
suggested  that  spatial  lagging  of  maximum  rj  and  i//  was  indicative  of 
horizontal  ’transitions’  of  the  scattering  layer  leading  to  changes  on 
depth,  amplitude  and  thickness  of  sub-surface  S  layers.  Based  on 
numerical  experiments.  Carpenter  et  al.  (2009)  proposed  that  drastic 
temporal  variations  of  skewness  and  kurtosis  can  reflect  dominance 
transitions  in  phytoplankton  communities. 

Detailed  analysis  of  lidar  profile  ’shape’  parameters  derived  from 
sub-surface  S  peaks  and  R^-ratio  statistics  showed  that  Rj  skewness 
was  sensitive  to  deepening  of  sub-surface  S  maximum.  In  general, 
lidar  profiles  characterized  by  a  single  S  sill  above  the  water  baseline 
had  more  left-skewed  R,  distributions  (i.e.,  more  negative  «//)  when 
this  shape  feature  was  located  far  from  the  sea  surface.  Based  on  light 
propagation  models,  Kerfoot  et  al.  (2005)  calculated  a  higher  R, 
when  Karenia  brevis ,  a  phototaxic  dinoflagellate,  was  concentrated 
near  the  ocean  surface  (0-2  m  depth).  This  R,  increase  was  likely 
reflecting  a  vertically  integrated  ocean  color  signal  dominated  by 
particles  smaller  than  K.  brevis.  As  large  cells  (>20pm)  of  K.  brevis 
migrate  down  during  morning  (5  to  12  am)  and  late  evening  (7  pm 
to  5  am),  R]  decreases  reflecting  a  re-distribution  of  the  surface  layer 
composed  by  ’large’  dinoflagellate-derived  particles  throughout  the 
water  column. 
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Fig.  A1.  Vertical  structure  of  density  and  optical  properties  during  August  2002.  a) 
Geographic  location  of  ship-based  sampling  stations  (asterisks),  b)  chi  o.  KPAr.  and 
seawater  density. 


In  NGOA  waters  during  summer  of  2002,  ocean  color  statistics 
based  on  Rn  ratios  were  not  better  predictors  of  g2  than  g5  and  g3  (see 
Eq.  (10),  Table  4).  This  limitation  was  somehow  evidenced  by  those 
studies  modeling  deep  chlorophyll  a  maximum  (DCM)  based  on 
satellite-derived  chlorophyll  a  concentration  (Millan  Nunez  et  al.. 
1997,  Sathyendranath  et  al.,  1995). 

As  expected,  parallel  changes  on  higher  statistical  modes  of  Rn 
ratios  and  £ 5  between  lidar  profiles  with  contrasting  vertical  shapes 
were  also  accompanied  by  modifications  on  lidar  k and  cv.  In  average, 
lidar  waveforms  with  sub-surface  Gaussian  peaks  had  circa  of  one 
order  of  magnitude  higher  power  spectrum  slope  of  5  (arithmetic 
average  of  k= 4.1 5  ±7.9,  (1  standard  deviation))  and  coefficient  of 
variation  per  bin  (250  m  horizontal  x  18  m  vertical)  with  respect  to 
those  lidar  profiles  with  monotonic  decrease  of  S  with  depth  (Fig.  3). 


Table  Al 

Correlation  (p,)  between  ‘shape*  parameters  describing  lidar  waveforms  with  a  single 
sub-surface  backscattering  maximum. 


Comparison 

ft 

n  * 

<N 

1 

C 

V 

P 

g2  versus  g, 

-0.48 

23 

—  232 

0.02 

%2  versus  g3 

0.28 

23 

135 

0.19 

g,  versus  g3 

035 

23 

1.72 

010 

J  n  is  equal  to  the  number  of  observations.  r(n  —  2)  is  the  5tudent  t-test  statistics,  p,  is 
significant  at  95%  confidence  level  when  P  0.05. 


Considering  all  lidar  profiles  obtained  during  August  17  2002  survey, 
our  arithmetic  average  k  was  0.61  ±0.44  (1  standard  deviation)  and 
significantly  lower  with  respect  to  the  mean  value  (1.49  ±0.03,  near 
surface  returns  neglected)  reported  by  Churnside  and  Wilson  (2006) 
based  on  FLOE  measurements  off  the  coast  of  Oregon  and  Washington 
during  July  of  2003.  As  pointed  out  by  Churnside  and  Wilson  (2006) 
most  of  our  h  values  indicate  that  lidar  backscattering  patches  were 
smaller  with  respect  to  those  patches  generated  by  turbulence  within 
the  inertial  subrange  or  Brownian  motion  (k=  1.67  to  2.0)  but  larger 
with  respect  to  spatial  fluctuations  equally  distributed  in  terms  of 
amplitude  (i.e.,  k=0  for  white  noise).  Lower  S  patchiness  in  this 
study  (i.e.,  higher  k)  with  respect  to  that  observed  by  Churnside 
and  Wilson  (2006)  can  be  attributed  to  a  multiplicity  of  factors  in¬ 
cluding  differences  of  plankton  communities,  water  stratification, 
primary  productivity,  and  vertical  migration  patterns  between 
sampling  locations. 


5.  Conclusions 

Relationships  between  vertical  structure  of  optical  properties  and 
surface  ocean  color  information  were  intensively  studied  during  the 
1990s  using  ocean  color  imagery  and  ship-based  measurements 
(Millan-Nunez  et  al.,  1997;  Sathyendranath  et  al.,  1995).  More 
recently,  the  optical  structure  of  the  ocean  intenor  has  been  mainly 
investigated  using  in-water  unmanned  sensors  (e.g.,  gliders.  Hodges  8? 
Frantantoni,  2009,  AUVs,  Sullivan  ct  al..  2009). 

In  this  work,  we  propose  an  alternative  approach  using  high 
resolution  active  (lidar)  and  passive  (ocean  color  radiometer)  optical 
remotely  sensed  data.  Briefly,  we  propose  the  use  of  low  and  high 
statistical  moments  of  Ri  to  detect  S  changes  along  the  vertical 
component.  Locations  characterized  by  maximum  values  of  S  toward 
the  sea  surface  were  generally  associated  with  Rj  spatial  distributions 
characterized  by  relatively  high  arithmetic  average,  standard  devia¬ 
tion,  kurtosis  and  skewness.  Conversely,  lidar  waveforms  with 
additional  sub-surface  backscattering  features  tended  to  have  lower 
arithmetic  average,  standard  deviation,  kurtosis  and  skewness  of  Rt . 
Whether  these  optical  functionalities  between  active  and  passive 
measurements  hold  in  other  marine  regions  or  other  periods  of  the 
year  is  a  matter  that  deserves  more  study 

In  this  study,  we  demonstrated  that  passive  and  active  optical 
measurements  can  be  connected  based  on  a  single  wavelength  (i.e,, 
A  =  532  nm)  even  though  the  apparent  lack  of  correlation  between  a  and 
MicroSAS-derived  optical  properties.  The  observed  relationships  between 
Rn  ratios  and  lidar  parameters  were  influenced  by  additional  variables 
that  would  require  additional  experiments  involving  characterization  of 
optical  targets,  different  lidar  polarizations  and  sunlight  geometries,  and 
simultaneous  in  situ  measurements  of  vertical  distribution  of  lOPs. 

Observed  optical  relationships  derived  from  this  study  represent  a 
preliminary  step  to  better  describe  thin  layers  and  DCM  dynamics, 
and  help  improve  QC  of  ocean  color  products,  at  least  in  coastal  waters 
such  as  those  observed  in  this  study,  waters  characterized  by  stratified 
conditions,  well  defined  lidar  backscattering  layers  (i.e.,  magnitude 
of  5  maximum  comparable  to  near  surface  values),  and  relatively 
high  vertically  integrated  chi  a  (>150  mg  m  2)  within  the  euphotic 
zone. 
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Appendix  A.  Ship-based  profiles  of  oceanographic  variables 

Shipboard  profiles  of  downwelling  photosynthetically  available 
radiation  (PAR)  were  conducted  aboard  FV  Laura  of  Kodiak  Island  and 
during  August  16-18  of  2002.  The  summer  survey  consisted  in  26 
vertical  casts  encompassing  a  depth  range  between  50  and  198  m 
(Fig.  Ala).  Surface  (Ed(PAR,  04-))  and  underwater  (Ed(PAR,0-))  PAR 
measurements  were  collected  with  a  Biospherical  Instruments,  Inc., 
PRR  2600  profiling  reflectance  radiometer  system.  Raw  Ed(PAR,0-) 
profiles  were  quality  checked  by  removing  spikes  and  values 
increasing  with  depth  caused  by  temporal  variation  of  near  surface 
(e.g.,  wave  focusing)  and  above  surface  (e.g.  cloudiness)  environ¬ 
mental  factors.  A  posteriori ,  downcast  Ed(PAR,0— )  profiles  were 
smoothed  every  meter  based  on  arithmetic  averaging,  and  ver¬ 
tically  diffuse  attenuation  coefficient  of  PAR  (KPAR)  was  calculated 
assuming  an  exponential  decrease  (neperian  base)  of  PAR  with  depth 
(e.g.,  Kpar(z2)  =  -  log(Ed(PAR,  0-)2l  /Ed(PAR,0-)z2)/(z2  -  z, )). 

Temperature,  practical  salinity,  chlorophyll  a  fluorescence  and 
water  depth  were  recorded  with  a  Sea-Bird  Electronics,  Inc.,  model 
25-143  CTD.  Pre-cruise  CTD  calibrations  were  performed  at  Sea-Bird 
lab  (temperature  accuracy  ±  0.002  °C,  practical  salinity  accuracy  ± 
0.005).  Magnitude  of  chi  a  was  estimated  from  fluorescence  mea¬ 
surements  (Seapoint  fluorometer,  excitation  wavelength  =  470  nm. 
emission  =  685  nm,  sensitivity  =  0.02  mg  m  3,  and  temperature  sta¬ 
bility  of  <0.2%fC).  The  Seapoint  sensor  was  calibrated  in  the  lab  with 
a  Turner  Designs  Model  10-AU  fluorometer  and  using  the  non- 
acification  method  (Welschmeyer,  1994).  In  general  terms,  the 
median  of  CTD-derived  variables  in  shelf  waters  of  Afgonak/Kodiak 
Islands  evidenced  stratified  conditions  characterized  by  a  developed 
pycnocline  in  the  first  10  m  of  the  water  column  (i.e.,  vertical  sea¬ 
water  density  changes  >0.02  kg  m  4)  associated  with  drastic  increase 
(>100%)  of  Kpar  above  that  depth  and  larger  chi  a  toward  a  water 
depth  of  7  m  (up  to  8.2  mg  m~3)  (Fig.  Alb). 
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